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ABSTRACT 

In this paper we describe an adaptive softening length formalism for collisionless 
N— body and self-gravitating Smoothed Particle Hydrodynamics (SPH) calculations 
which conserves momentum and energy exactly. This means that spatially variable 
softening lengths can be used in N — body calculations without secular increases in en- 
ergy. The formalism requires the calculation of a small additional term to the gravita- 
tional force related to the gradient of the softening length. The extra term is similar in 
form to the usual SPH pressure force (although opposite in direction) and is therefore 
straightforward to implement in any SPH code at almost no extra cost. For N— body 
codes some additional cost is involved as the formalism requires the computation of 
the density via a summation over neighbouring particles using the smoothing kernel. 
The results of numerical tests demonstrate that, for homogeneous mass distributions, 
the use of adaptive softening lengths gives a softening which is always close to the 
'optimal' choice of fixed softening parameter, removing the need for fine-tuning. For 
a heterogeneous mass distribution (as may be found in any large scale N— body sim- 
ulation) we find that the errors on the least-dense component are lowered by an order 
of magnitude compared to the use of a fixed softening length tuned to the densest 
component. For SPH codes our method presents a natural and elegant choice of soft- 
ening formalism which makes a small improvement to both the force resolution and 
the total energy conservation at almost zero additional cost. 

Key words: hydrodynamics - methods: numerical - gravitation - methods:./V-body 
simulations 



1 INTRODUCTION 

This paper is concerned with the question of how best to rep- 
resent the gravitational force when simulating self gravitat- 
ing systems using particles. The simplest of such systems is 
a collection of stars which is usually replaced by a very much 
smaller number of computational particles. A more compli- 
cated example is the simulation of self gravitating gas with 
or wi thout a stellar co mponent, using the particle method 
SPH <Monaghadl2005h . 

Provided the number of computational particles is suf- 
ficient to resolve the important dynamical scales the simula- 
tion can give satisfactory results for most quantities though 
the slow relaxation of a galaxy is number dependent. In 
the case of A''— body simulations it is, however, necessary 
to soften or smooth the forces between pairs of particles so 
that the binary collisions of the computational particles will 
not cause numerical artifacts. 

The simple Plummer form of the softening where the 



force F(r) between a particle pair with masses m a and mt 
separated by distance r is 

F ^) = -G {r2 + h2) . A/2 , (1) 

and h is the softening length. lDermenl<|200ll) . amongst others 
fe.g. lDver fc IrJl993l) has shown that a better choice is to use 
Kernel smoothing with kernels W(r, h) that have compact 
support. The softened force at a due to b then takes the form 



F = —G Y'" / W{r,h)r 2 dr. (2) 
r Jo 

Provided h is constant Poison's equation shows that the soft- 
ening is equivalent to calculating the local gravitational force 
on a point particle a due to a density 

p(v a ) = m b W{r a - v b ,h), (3) 

or, when there is a collection of particles contributing to the 
force on particle a, the density is 
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p(r a ) = 22 mbW(r a -r b ,h) 



(4) 



which is identical to the SPH density estimate. 

Kernels with compact support are zero beyond some 
specified distance proportional to the length scale h, and 
the pair force then has the correct value for the two sets of 
real particles represented by two computational particles. 

The usual practice is to use a fixed value of h for all of 
the iV— body particles. A key issue which arises in this con- 
text ,jmd_^iclihasj3e£rit]^ 

ies iMerrittl |l996t iRomeol Il998t lAthanassoula et alJ l2000t 



lDehnenll200ll : lRodionov fc Sotnikovall2005l) . is the 'optimal' 
choice of softening length, for too small a softening length 
will result in noisy force estimates, whilst too large a value 
will systematically 'bias' the force in an unphysical manner. 
In general, however, the 'optimal' choice depends on partic- 
ular sy stem under investigation and may not b e known a 
priori iAthanassoula et alJl200oT ). iDehnenl 12001"! quantifies 
the errors arising from both Plummer and kernel softening of 
the above form. In all cases he finds the accuracy is improved 
if h is allowed to vary according to the local particle number 
density n in such a way that h is smaller when n is larger. 
Typically h tx 1/n 1 '' 3 . For self-gravitating SPH calculations 
using a softening length which diffe rs from the smoothin g 
length can lead to unphysical results feate fc Burkertll997l) . 

To retain conservation of linear and angular momen- 
tum it is necessary to use a symmetric form of F so that 
each particle in a pair interaction experiences an equal but 
opposite force. This can be achieved by using, for example 
h — ^(hi+hj) in place of h in (|2J. However, because the soft- 
ening length then varies in space the total energy of the sys- 
tem will not be conserved. The errors are often not large but 
can lead to secular increases in the total energy of the sys- 
tem, destroying the phase-space conse rvation which is cru- 
cial for accurate N —body simulatio ji_jH^riio t uisLfc Barnes! 
Il990l: IPehnenlEoOll : iRodionov fc Sotnikoval^OOftl . 

In this paper we show how a Lagrangian for a self grav- 
itating gas can be devised which has the softening of the 
force and the variation of h built in. The advantage of using 
a Lagrangian is that, provided it is constructed correctly, the 
conservation laws are automatically satisfied. In particular 
the conservation of energy and momentum is exact though, 
in practice, the accuracy is determined by the time stepping 
algorithm. The new equations of motion have an extra term 
in addition to the standard SPH and gravity terms. It is 
this term which guarantees energy conservation. We apply 
our algorithm to both static and dynamic problems. In some 
cases the new equations give results which are very similar 
to results obtained previously, but in some cases the results 
are very much improved. 



2 KERNEL SOFTENING 

A general for mulation for force softening was given by 
IDehnenl J200ll) and we use a similar formulation here. The 
modified gravitational potential per unit mass may be writ- 
ten in the form 



JV 



$(r) = -G^m b ct>(\r- 



■ n\,h) 



(5) 
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Figure 1. The functional form of the modified potential (-), grav- 
itational force and the density profile using the cubic spline kernel. 
For r/h>2 the smoothing is zero and the potential and force are 
exact. 



where (j> is the softening kernel which is a function of the 
particle separation and the softening length h (we use h to 
denote the softening length since it corresponds with the 
smoothing length used in the SPH density estimate). The 
kernel determines the functional form of the modified grav- 
itational force law. For example, in the case of Plummer 
softening the softening kernel is given by 



(r,h) = - h 



1 + 



(s)' 



-1/2 



(6) 



Neglecting the spatial variation of h the force estimate 
corresponding to © is given by 



F(r) = -V* = -G^m50'(|r-rfc|,ft) 



r - r b 
|r - rJ 



(7) 



where <f>' — dcj)/d\v — | . The underlying smooth density 
field can be obtained from Poisson's equation 



V 2 $ = 4ttGp, 



giving 



p(r) = J2 m >>W(\r-r b \,h) : 



(8) 



(9) 



where the density kernel is related to the softening kernel 
according to 



W(r) 



d 



(10) 



4nr 2 dr 

The kernel density given by © corresponds to the mass dis- 
tribution of each particle being smoothed. Readers familiar 
with SPH will notice that JUJ corresponds to the density es- 
timate used in SPH calculations, where W is the usual SPH 
smoothing kernel. 

In general, the functional form of the softening kernel 
may be specified for either the potential term <f>, the force 
evaluation term <f>' or W. In each case the corresponding 
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kernel for the other cases may be determined by a straight- 
forward integration or differentiation. For example, in N- 
body codes, it may be preferable to choose a kernel pri- 
marily for the force evaluation, from which the functional 
form of the potential and density kernel can be derived. In 
SPH the kernel is primarily used for the density estimate, 
where the most commonly us ed form is the cubic spline of 
iMonaghan fe Lattanzlcl il985ft 



W(r,h) = — l J(2- 9 )«, 



< q< 1; 

1 < q < 2; 
q > 2. 



(11) 



where q = r/h. The corresponding force kernel is given by 

,/ 47T 



Wrdr', 



(12) 



the functional form of which is given for the cubic spline 
in Appendix^ The softening kernel for the gravitational 
potential may be calculated from the force kernel using 



= / Fdr, 



(13) 



the form of which is also given in AppendixEJfor the cubic 
spline. For general kernels I13H combined with 1121 can be 
integrated by parts to give 



(r, ft) = An 



Wr' 2 dr' + / Wr'dr' - 



Wr'dr' 



where the last term is the constant of integration, deter- 
mined by the requirement that <fi — > as r — > oo (note that 
we have also assumed a kernel with compact support of size 
2ft). 

The modified potential, force functions and the density 
kernel are shown in Figure^for the cubic spline. The reader 
should note that whilst we use the cubic spline as an example 
throughout this paper, the results derived in the following 
sections are quite general and any smooth ing kernel may be 
used (including any of those suggested by lDehnenll200lf) . 



3 LAGRANGIAN FORMULATION 

The Lagrangian describing a self-gravitating gas is given by 



1 

L = m f> ( g v b ~ ~ "'• 

where $ is the gravitational potential @ and u is the ther- 
mal energy per unit mass. The equations of motion may be 
obtained via the Euler-Lagrange equations 

d ( 8L\ dL „ 



giving 
dv a 



dL 

dv a 



(16) 



The advantage of using a Lagrangian to derive the equa- 
tions of motion is that, provided the Lagrangian is sym- 
metrised appropriately, momentum and energy conservation 
are guaranteed. Variational principles have been used exten- 
sively to derive conservative SPH formalisms for relativistic 



fluid dyna mics jMonaghan fc Price ! Eool . magnetohydro- 



dynamics dPrice fc^^onarfianf 



20041) and in the case of a 



spatially variable smo othing length l|SDringel fc Hernauisd 
l2002l lMonagharJl2002l) . 

An adaptive softening length formalism may be derived 
by writing the gravitational part of the Lagrangian in the 
form 



b 
G 



2 EE 

b c 



mbm c (f>bc(hb) , 



(17) 



where <j>b c refers to 0(|rf, — r e [). Swapping indices in the 
double summation shows that 1171 is equivalent to averaging 
the softening kernels in the form 



G 



EE 



mi,m c 



4>bc{h b ) + 4>bc{h c 



(18) 



The derivative of l)17|l is given by 



dL a 



dr a 



where 

d\r bc \ 
<9r a 



d(j) bc {h h ) 



+ 



d\r bc \ 
d(j>bc(hb) 



dh b 



o\r bc \ 
<9r a 

dh b 
dv a 



r b - r c 
\r b - r c \ 



(5 ba — Sea)- 



(19) 



(20) 



We relate the smoothing length to the particle co-ordinates 
assuming h = h(p), using 

dh b _ dh b dp b 
dr a dp b dv a ' 

where p is the density calculated by a summation over neigh- 
bouring particles in the form 



(21) 



pa = m b W (\r a - r b \,h a ) , 



(22) 



where W is the density kernel. The relationship between h 
and p means that this is a non-linear equation for both h a 
and p a which can be solved self-consistently for each particle. 
The iterative method we use for doing so is described in 
detail in iJ4.2l The spatial derivative of 1221 is 



(14) £a = _l 

dr a 



1 dW b d(hb) r , 

TT- > m, d r (O ba — dda) ■ 

Q b ^— ' dr a 

d 



(23) 



where W is the density kernel and SI is a term accounting 
for the gradient of the smoothing length given by 



dh a ST^ dWab(ha) 

1_ a — / mb HI 

op a ^ ah a 

b 



(24) 



have 
dL„, 



Using lEH and ll^3Tl in iTHJIl and simplifying, we 

4>'ab{ h a) + (Pab{h b 



dv a 



E 

6 



-m a m b 



^E 



m b - 



( a dW ab (h a ) 



2 \ Sl a dr a 



+ 



r a - r b 

|r a - r b \ 

(b dWab(hb) 



(25) 



dv a 



The quantity £ is defined as 
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dhc 

dpa 



Ed(j) ab (h a 
m b 

6 



dh a 



(26) 



where d(j>/dh can be tabulated (or calculated) directly for 
the particular smoothing kernel used. For the cubic spline 
the expression is given in Appendix lAl 

The derivation of the SPH pressure force from the ther- 
mal energy term in the Lagrangian 11411 in the case of a spa- 
tially variable smoothing length has been described in detail 
elsewhere (e.g.lSpringel fc Hernauisi2002tlMonaghadl20o3: 
IPrice fc MonaghanB2004l) and we simply use the result here. 
The final equations of motion take the form 



dt 



-G22 m b 

b 



<t>'ab(ha) + <j>' ab (h b ) 



Gy 

b 



b 



nib 



(a dW ab (h a ) 

£l a dr a 

Pa dW ab (h a ) 

plQ a dr a 



+ 



r a - r b \ 

(b dWabjhb) 
fl& dv a 

Pb dWg b (h b 



(27) 



dr a 



The first term in l|270 corresponds to the softened gravita- 
tional force. The second term is present only in the case 
of adaptive softening lengths and it is the incorporation of 
this term which restores the energy conservation. The third 
term is the usual SPH pressure force allowing for a spatially 
variable smoothing length. The terms fi and £ required in 
the adaptive softening term (and for also in the pressure 
force) are easily calculated alongside the density summation. 

The additional adaptive softening length term can be 
seen to have the same form as the pressure force, with the 
quantity P/p 2 replaced by £. Notice however that £ is, for 
positive kernels, a negative definite quantity and therefore 
that the adaptive softening term acts in opposition to the 
usual pressure term (that is, in the direction of increasing the 
gravitational force). This is in line with the recent findings of 
iHubber et alJ J2006) ; that SPH always underestimates the 
gravitational force at low resolution. In fact they suggest 
adding an additional contribution to the gravitational force 
based on the 'self-gravity' of an SPH particle. The new term 
derived above provides a similar contribution without the 
need for ad-hoc prescriptions. 

Alternative formulations of the adaptive softening for- 
malism given above are possible by symmetrising the La- 
grangian in different ways. We use the formulation given 
above since it is simple and efficient to implement. As an 
example, we derive an alternative version based on the aver- 
age softening length in Appendix^] Whilst the force derived 
using the average softening length is very similar to 12711 but 
in keeping with the average h used in th e variable smooth- 
ing length SPH formalism of lBend (ll990T) . it has the practi- 
cal disadvantage that we do not use the average smoothing 
length elsewhere in the calculations, neither in the density 
summation (since the density for particle a is calculated us- 
ing only ha) nor in the SPH force term (see 1271 . Thus the 
calculation of the quantity £ IB6> for each particle requires 
the calculation of the kernel not only using h a (for the den- 
sity and fi) but additionally using h a b which is not only 
inefficient but also rather inelegant in the numerical code. 
Thus, we do not use the average h formalism in this paper. 

A further possibility, not examined in this paper, would 
be to use the product of the softening kernels in the force 



evaluation. The situation is complicated slightly in this case 
as the product form must either be used in the potential or 
force but not both. In any case the differences in the force 
evaluated using different symme trised form s are v ery small. 
The suggestion put forward bv lDver fc id (ll993T> that the 
force should be symmetrised by considering two overlapping 
spheres may also be used in a similar manner to derive an 
energy conserving formalism, but it is not cle ar that there i s 
any advantage to be gained by doing so (see Dehnen 2001). 

For reference the consistent forms of the continuity and 
internal energy equations for SPH simulations are given by 



dpa 1 ( \ 

-dT = ir a 2^ mb( - Va - Vb) 



dt 

and 

du a 
~dt 



Pa 

Sl a pl 2^ 
b 



dWg b (h a ) 

dv a 



m b (v a - v b ) 



dWg b (h a ) 
dv a 



(28) 



(29) 



where v is the particle velocity. The continuity equation can 
be used to make a starting guess for the h iteration pro- 
cedure used to determine the density (described in the fol- 
lowing section). An alternative to using the internal energy 
equat ion is to evolve the entrop y as an independent vari- 
able llSpringel fc Hernau ist 2002) which is possible for ideal 
equations of state. 



4 NUMERICAL TESTS 

We test the adaptive softening length formalism derived 
in the previous section using three examples . The first 
(MJil is a series of static tests used by iDehnenl (120011 ) and 
Atha nassoula et alJ a2000r) in order to estimate the force er- 
rors associated with softening formulations. We also consider 
a dynamic version of one of these tests in order to study the 
energy conservation properties of our new method ( jEQfl . 
The second example ( H4.4H involves self-gravitating SPH and 
the static structure and dynamical oscillation of a polytrope. 



4.1 Errors 

In the static halo tests, we calculate the Average Square 
Error (ASE) in the gravitational force according to 



ASE = ^J2\fi ~ fe X act{xi)\\ 



(30) 



where /» is the force on particle i, N is the particle number 
and C is a normalisation constant. Unless otherwise specified 
we use C — 1/ where f m ax is the maximum value of the 
exact solution. The Mean Average Square Error (MASE) is 
then the mean over all realisations, 



/ N 

MASE = — / ^ l/< - /«««*(* 
\ i=i 



(31) 



We choose this quantity rathe r than the Me an In tegrated 
Squar e Error (MISE) used by iMerritd dl996ft and IDehnenl 
feOOll) . given by 



MISE - 



C_ 

M 



P( x )l/( x ) - /«Mct(x)| dx 



(32) 
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where p(x) is the true density at a point x, /(x) is the force 
calculated at that point from the N— body distribution and 
M is the total mass. Calculation of the MISE is complicated 
by the need to integrate along radial grid points (involving 
calculation of the force at positions ot her than particle po- 
sitions) and lAthanassoula et al.l i200(J) find little difference 
between their results using MASE or MISE error measures. 
In our case the correction terms derived in depend on 
a particle's own density estimate, so it makes sense to cal- 
culate errors only at particle positions (that is, using the 
MASE estimate). 

The reader should bear in mind that, using either the 
MASE or MISE as defined above, the total error tends to 
be dominated by the regions containing the largest forces. 
This can be somewhat misleading in comparing adaptive 
softening with fixed softening, as the fixed softening length is 
generally chosen to minimise the error in the densest regions, 
where the adaptive softening will not show a large difference. 
An example is given in >I4.3.2I where a two-halo system is 
setup and we explicitly show the contribution to the MASE 
from each halo, even though the total error is dominated by 
the densest component. 

4.2 Setting the softening length 

The method we use for setting th e softening length 
is identical to the met hod used by iPricd (12004?) (see 
iPrice fc Monaghanlliooi) for setting the smoothing length 
in SPH calculations. A sim ilar method is also used by 
ISpringel fc Hernauistl i2002T) and hence also in the publicly 
availa ble GADGET-2 code for N-body and SPH iSpringell 
120051) . The idea is to regard the smoothing length as a func- 
tion of density via the relation 



u -1/3 
ha OC p a 



or more specifically 

1/3 

, / m a 

h a =ri\ 

V Pa 



(33) 



(34) 



where m is the particle mass, p is the mass density and 
r\ is a dimensionless parameter which specifies the size of 
the smoothing length in terms of the a verage particle spac- 
ing (similar to the parameter e used bv iDehnenl EoO 1) . The 
derivative is given by 



(35) 



dh a _ _ h a 
dp a Spa ' 

An equivalent interpretation of l|34|l is that the mass 
contained within a smoot hing sphere is held constant 
jSnringel fc HernmiisdEool. that is 

4 3 

--ir(aha) pa = const = m a N nE i g h, (36) 

where a is the compact support radius of the kernel (= 2 for 
the cubic spline) and N„ e i g h = ^n(arj) 3 may be used as an 
approximate measure of the number of neighbours contained 
within a smoothing sphere. Unless otherwise specified we use 
r\ = 1.2 in the variable smoothing/softening length formula- 
tions used throughout this paper, which in three dimensions 
is equivalent to ~ 60 neighbours. 

For a pure N-body simulation using unequal mass bod- 
ies, it may be advantageous to use a number density rather 



than the mass density for setting the softening length. The 
resulting gravitational force in that case is identical to the 
first two terms in 1271 with mass density replaced by number 
density. In this paper we assume a mass density dependence 
consistent with SPH simulations. 

The density is calculated by a direct summation over 
the particles in the form 122H which, via the relation 1341 
be comes a non-lin ear equation to be solved for both h and 
p. iDehnenl i200ll) suggests that even a rough approxima- 
tion of the (number) density is sufficient for the purpose 
of adapting the softening length in A^— body calculations. 
A similar argument may be made for setting the smooth- 
ing length in SPH calculations. In both cases, however, the 
situation changes once the gradient terms are incorporated 
into the equations of motion (as in this paper) since these 
terms are calculated on the basis of the h(p) (or h(n)) re- 
lation and may therefore introduce substantial inaccuracies 
into the solution if the density and smoothing (or softening) 
length are far from being consistent with 1348 . 

From a practical point of view obtaining a self- 
consistent solution to (1341 and 1221 is a relatively straight- 
forward root-finding problem. The function to be solved may 
be written in terms of either h or p. Written in terms of h 
we have 



f(h) = 0, 
in the form 

pa{h a ) — Psum{h a ) = 0, 



(37) 



(38) 

where p a is the density consistent with the current smooth- 
ing length h a calculated from the relation I34H and p S um is 
the density calculated using h a from the summation over 
neighbouring particles 12211 . We use a Newton- Raphson it- 
eration method, ie. 



ha 



= h a 



f(ha 



f'(ha) 

where the derivative of l|380 is given by 



f'(ha) = 



d£a_ ST^ dWgbjhq) _ 3p a 

dh a ^ dh a h a 

b 



n . 



(39) 



(40) 



We find this method to be efficient and cost effective, par- 
ticularly since the quantity Q (defined in i'24l ) is already 
calculated alongside the density summation for use in the 
equations of motion. 

Convergence is determined for each particle individu- 
ally according to the criterion \h„ ew — h\/ho < e where ho 
is the smoothing (or softening) length at the start of the 
iteration procedure and typically we use e = 10~ 3 . We find 
that it is more efficient to perform the iterations by loop- 
ing over the particles as the outer loop and iterating each 
particle at a time to convergence. We also find that it is no 
longer efficient to store a global neighbour list for all par- 
ticles but rather to perform a neighbour search on-the-fly 
(e.g. using a treecode), recalculating where necessary and 
being stored only for the particle being iterated. This also 
represents a significant reduction in memory requirements 
for SPH calculations. 

The Newton-Raphson iterations work extremely well 
provided that the initial estimate of h is reasonably close 
to the actual solution. This is almost always the case in the 
calculations since there is only a small change in p between 
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timcsteps. The only problems which may arise are in the 
first iterations on the initial conditions where h and p may 
be far from the relation l)34jl . For this reason it is useful to 
revert to a bisection scheme (which is guaranteed to con- 
verge) in the case where the Newton-Raphson iterations do 
not converge (we set the limit for this as > 20 iterations). 

In terms of cost the density iterations add only a small 
amount of extra work to SPH calculations. The exact work 
required depends on the nature of the simulation since more 
iterations are required when the density is changing rapidly. 
However the scheme is very efficient since iterations are only 
performed on the subset of particles whose densities are 
changing. The scheme can be made still more efficient by 
predicting an initial guess for h in the time evolution scheme 
using the time derivative 



(41) 



dh dh dp 
dt dp dt ' 

where for the summation 1221 the time derivative is given 
by (J2HJ- Using a prediction step we find that in general (al- 
though dependent on the dynamics of a particular simu- 
lation) only a small fraction of the particles require extra 
density calculations and that these particles then converge 
rapidly (in ~ 2 — 3 iterations) . 



4.3 iV-body tests 

4-3.1 Isolated haloes 

The first test we perform is to compare the (softened) grav- 
itational force to the exact force given an analytic density 
profile (corresponding to typical structures formed in cos- 
mological simulati ons or used as i nitial conditions in galaxy 
models). Following iDehnenl (|200ltl we consider two different 
density profiles - Plummer spheres and Hernquist models. 
The density profile for the Plummer spheres are given by 



p(r) 



3GMr 2 



47r(rl +r 2 ) 5 /2' 



(42) 



where M is the mass and r s is a parameter determining the 
concentration of the halo. The corresponding gravitational 
potential and force are given by 



$(r) = 
and 

$'(r) = 



GM 



( r 2 +r 2)l/2> 



GM 

( r 2 +r 2)3/2 



(43) 



(44) 



The cumulative mass profile for the Plummer sphere is given 
by 



M(r) 



GMr 



(45) 



( r ,2 + J .2)3/2- 

The lHernauisd(ll99d) model consists of a density profile 
given by 

GMr B 



p(r) 



27rr(r s + r) 3 



(46) 



The density and gravitational forces in this model are more 
difficult to resolve since the density profile is cusped near 



the origin. The potential and force are given by 
GM 



$(r) = - 
$'(r) 



(r s + r) ' 
GM 



{r s + r) 2 



whilst the mass profile is given by 



M(r) = 



GMr 2 
(r s + r) 2 



(47) 
(48) 

(49) 



The density profiles in each case are set up in the 
usual manner choosing three random deviates (xi, 22, £3) 
uniformly on (0, 1). The first is used as a position in the 
mass profile from which the radial coordinate is determined 
by rearranging M(r) to give r(m) where m is the mass frac- 
tion. In practice we use only mass fractions < 0.99 in order 
to prevent isolated particles being placed at extremely low 
densities. The second random number xi is used to give a 
random azimuthal angle ip — n(2x2 — 1) whilst the third 
is used to give a spherical angle 9 via the transformation 
9 — cos -1 (2x3 — 1) (necessary to prevent the distribution 
from clumping towards the poles). The result is a particle 
distribution which closely mirrors the analytic density pro- 
file although with errors decreasing like 1/y/N due to the 
Monte Carlo nature of the distribution. 

In the numerical simulations we use units of mass 
[M] = 1, length [R] = 1 and time [r] = (GM/R 3 )~ 1/2 . 
In these units GM — 1 such that the gravitational constant 
does not appear in the numerical equations. Correspond- 
ingly, force and energy (both per unit mass) are measured 
in units of GM/ R 2 and GM/R respectively. For calculation 
of the mean error we compute 3 x 10 6 /N realisations for a 
halo of N particles. 

The MASE calculated for Plummer haloes of N = 10 2 , 
10 3 , 10 4 and 10 5 particles with M = 1 and r 3 = 1 are 
shown in Figure |2] The left panel shows the results using 
a fixed softening length, comparing both Plummer (solid 
lines) and cubic spline (dashed lines) softening kernels. The 
right panel shows the results using (cubic spline) adaptive 
softening, with (dashed line) and without (solid line) the 
energy-conserving term, also showing the variation with the 
"Number of Neighbours" , by which we mean the parameter 
N ne igh defined in ll3lH . 

Some general features are worth pointing out. Firstly, 
using a fixed softening length, there is a large variation in the 
total error depending on the choice of softening length (left 
panel). For softening lengths too small, the error is domi- 
nated by noise, reaching a maximum value once the soften- 
ing length is smaller than the smallest particle separation. 
For softening lengths too large, the error is dominated by 
the bias in the force introduced by the softening procedure. 
For some intermediate choice of softening length, there is 
a balance between noise and bias which produces a mini- 
mum error. This giv es rise to the c oncept of 'optimal soft- 
en ing' introduced bvlMerritd ( 1996) an d examine d in detail 
bv lAthanassoula et alJ feOOfl) and fPehnenl (1200 if) , whereby 
the softening length can be 'fine tuned' for a particular sim- 
ulation for minimum error. In principle this means that, for 
every iV— body calculation, there is an optimal choice of soft- 
ening length. The problem, demonstrated in the left panel 
of Figure [5] is that this 'optimal' choice not only depends 
on the parameters of the problem (for example, Figure [5] 
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Figure 2. MASE calculated for 3 X W^/N realisations of an 
isolated Plummer sphere with M = 1, r s = 1 and N = 10 2 , 10 3 , 
10 4 and 10 s particles. The left panel shows results using a fixed 
softening length, comparing Plummer softening (solid line) with 
cubic spline softening (dashed line). The right panel shows results 
using adaptive softening with (dashed line) and without (solid 
line) the new energy-conserving term. To guide the reader our 
fiducial choice of N n ei g h ~ 60 is indicated by the open circles. 



Hernquist model 



fixed plummer 
fixed cubic spline 



adaptive cubic spline 

adaptive + extra term 




0.01 0.1 1 

softening length 



50 100 200 500 
Nneigh 



Figure 3. As in Figure [5] but for an isolated Hernquist model 
with M = 1 and r s = 1. 



shows that the optimal ch oice changes with resolution, but 
lAthanassoula et alJ ll2000l) also discuss the dependence on 
the shape and degree of central concentration of the halo), 
but may also be different for different components of the 
same simulation. 

By contrast, use of adaptive softening (right panel) 
shows only a weak dependence on the choice of the (adap- 
tive) softening parameter (provided that the neighbour num- 
ber is small compared to the total number of particles). To 
guide the reader, our fiducial choice of Nneigh — 60 (given 
by 77 = 1.2 in (1341) ^1 is marked by the open circle in each case. 
Making this reasonable choice in all cases gives a softening 
which (according to the MASE estimate) is close to the op- 
timal choice of fixed softening. The exception is perhaps for 
the 10 4 particle halo, where the optimal choice of fixed soft- 
ening gives a MASE ~ 2.5 x lower than that using adaptive 
softening with energy conservation. However, changing the 



fixed softening length up or down by a factor of two in either 
direction (ie. not a large range if the optimal choice is not 
known a priori) means that even in this case that adaptive 
softening, even with the energy conservation term added, 



Comparison of the dashed (Plummer kernel) and solid 
lines (cubic spline kernel) in the left panel of Figure con- 
firms the conclusion reached by several authors that use of 
the cubic spline kernel is advantageous over the standard 
Plummer softening. In particular the optimal error for each 
halo is reached at a higher softening length for the cubic 
spline kernel and the slope in the error curve at high soften- 
ing lengths is steepened, demonstrating that the cubic spline 
reduces the bias in the force estimate. This is a result of the 
compact support of the cubic spline kernel, which gives a 
force with zero bias outside of the kernel radius (see discus- 
sion in 0. 

The right panel of Figure [5] also shows the influence of 
the new energy-conserving term in 1271 on the force errors 
in a static configuration. This term appears to increase the 
noise but also lower the bias in the (adaptively softened) 
force estimate, meaning that the total error is greater for 
smaller neighbour numbers but lower for larger neighbour 
numbers. We attribute this to the fact that the extra term 
is related to the gradients in softening length: Where these 
gradients are spurious (due to noise), the extra term may 
increase the total error. Where the gradients are due to ac- 
tual gradients in the density, the extra term correspondingly 
leads to a more accurate force estimate. This conclusion is 
also borne out by the results using Hernquist models (right 
panel of Figure |3J . Here the extra term leads to a smaller 
MASE (compare the dashed and solid lines in the right hand 
panel of Figure [3J at a lower N ne i g h than in the Plummer 
case (dashed vs solid lines in the right hand panel of Fig- 
ure|2J . The density profile in the Hernquist model is strongly 
cusped near the origin, meaning that any improvement in 
the resolution of density gradients (e.g. from the new term) 
tends to improve the error estimate. 

The Hernquist model was computed using M = 1, r s = 
1 and N = 10 2 , 10 3 , 10 4 and 10 5 particles. The results using 
fixed Plummer (solid lines) and cubic spline (dashed lines) 
softening on the Hernquist model are shown in the left hand 
panel of Figure [3] The differences between Plummer and 
cubic spline softening are much smaller in this case than for 
the Plummer spheres (Figure apart from a factor of ~ 2 
difference in the optimal choice of softening length for each 
kernel (that is, the optimal softening length for Plummer 
softening is approximately ~ 1/2 of the optimal value using 
cubic spline softening). 

These tests demonstrate that, for an isolated halo, the 
use of adaptive softening gives force errors which are close 
to optimal. Whilst there is not a significant improvement in 
the MASE compared to the use of an optimally-chosen fixed 
softening length, the use of adaptive softening removes the 
need for such fine-tuning. Furthermore it may not be possi- 
ble to find a softening which is 'optimal' for all components 
of a simulation. In the following section we consider such an 
example, where the use of adaptive softening shows a clear 
improvement. 
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4-3.2 Two Plummer spheres 

Next we consider two Plummer spheres placed at a fixed 
distance from each other, of equal mass but where one halo 
is much denser than the other. This situation may be rep- 
resentative of two haloes present in a typical cosmological 
N— body simulation or in a simulation of galaxy dynamics 
where more t han one galaxy is present. A similar test was 
considered bv lAthanassoula et al.l ll200(J) where a variety of 
mass ratios were also examined. Here we simply choose one 
representative case. 

Both haloes are Plummer spheres, setup as described in 
H4.3.1I We use equal mass spheres with M — 0.5. The first 
sphere is placed at the origin, with concentration parame- 
ter r 3 — 1 whilst a second sphere with r a =0.1 (ie. much 
denser) is placed some distance away at [as, y, z] = [10,0,0]. 
The MASE is calculated un-normalised in this case, that 
is with the normalisation factor C — 1 in order to make a 
meaningful comparison between the errors in each compo- 
nent. 

The MASE resulting from 300 realisations of this con- 
figuration using a total of 10,000 particles (5000 per sphere) 
is shown in Figure 2] (solid line), using fixed cubic spline 
softening (left panel) and adaptive softening with the en- 
ergy conservation term included (right panel), showing the 
variation with softening length in the former case and N nB i g h 
in the latter (where again the open circles correspond to our 
fiducial choice of N„ eigh ~ 60). The total MASE (solid line) 
is completely dominated by the densest component, with re- 
sults comparable to those shown in Figure [5] However, we 
also plot the contribution to the total MASE from the less- 
dense component (that is, the r s = 1 sphere at the origin) 
(dashed line). 

The problem with the use of a fixed softening length 
in a general N— body simulation is evident from Figure |1J 
namely that the 'optimal' choice of softening length differs 
for each component. The choice which minimises the errors 
in the densest component produces errors in the least dense 
component that are over l| orders of magnitude larger than 
the optimal choice of softening for that component. Con- 
versely choosing the softening which is optimal for the least 
dense component produces a bias in the force estimate in the 
densest component leading to a MASE ~ 2 orders of mag- 
nitude larger than the optimal choice of softening length in 
the dense component. Usual practise is therefore to choose 
the softening which minimises the softening in the densest 
component (s) (since this represents the largest contribution 
to the total error). However frequently one is interested in 
the properties of both (or all) components in an N— body 
simulation. This leads naturally to a need for adaptive force 
softening. 

The results using our adaptive softening length formal- 
ism (including the energy conservation term) are shown in 
the right hand panel of Figure 0] For both components the 
resulting error is, as previously, close to the optimal choice 
of fixed softening, but here the softening is close to optimal 
for both components. This means that the force errors in the 
least dense component are ~ 1 order of magnitude smaller 
than using a fixed softening length tuned to the densest 
component. 
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Figure 4. Mean averaged square errors in the gravitational force 
calculated for 300 realisations of a configuration involving two 
Plummer spheres. The total MASE is given by the solid line 
whilst the contribution to the total MASE from the least dense 
component is given by the dashed line. Results using a fixed cu- 
bic spline softening, varying the softening length, are shown in 
the left panel. The right panel shows the results using our adap- 
tive softening formalism (including the energy conservation term), 
varying the N neig h parameter. The open circles correspond to our 



fiducial choice of TV, 



neigh 



: 60. 



4-3.3 Halo relaxation 

An extension to the static halo test is to examine the dy- 
namic influence of the energy-conserving term. The ini- 
tial conditions for this test are a Plummer sphere with an 
isotropic velocity distribution corresponding to a (dynamic) 
steady state. The particle velocities are setup i n the manner 
described bv lAarseth. Henon fc Wielenl lll974l) : The distri- 
bution function is 



/(r,v,0) 



24y/2 
7vr 3 ( 





-E) 7 ' 2 E < 0, 
E > 0, 



(50) 



where /(r, v, t)drdv is the total mass of particles with posi- 
tion r and velocity v at time t and E is the energy per unit 
mass of a body: 



E 



1 2 V ' 2+<S> - 



(51) 



The distribution function is sampled by scaling the velocities 
in terms of the maximum velocity at r, ie. the escape velocity 



= v^2$ = 



2GM 



( r 2 +r 2)l/4 



(52) 



Writing q = v/v esc , from (1 -~> 1 p and l|50[l the probability dis- 
tribution for q is proportional to 



(53) 



where \q\ < 1. This distributi on is sampled usin g the Von 
Neumann rejection technique l|Press et alJ ll992*l: two uni- 
form random deviates X4 and x 5 are drawn. Noting that 
g(q) is always less than 0.1 (since \q\ < 1) we adopt q — X4 
if O.IX5 < g(q), otherwise a new pair of random numbers is 
tried until the inequality is satisfied. The velocity modulus v 
is obtained using 152H and, using two more uniform random 
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deviates xq and Xr, the velocities are given by 
= (l-2x 6 )v, 
= 

= \ru 2 



v x 

Vy 

Vz 



- Vx cos (2-11x7), 



- Vx sin (2tyx7 



(54) 



The halo is evolved forwards in time using a standard 
second order leapfrog integrator with a global timestep con- 
trolled by the condition 



At = 0.15 



1/2 



(55) 



where h is the softening length, / is the force per unit mass 
and the minimum over all particles is used. 

The energy conservation during the evolution of the 
equilibrium halo is shown in Figure |5] Using adaptive soft- 
ening without the additional term (solid line), fluctuations 
in the energy are observed from the changes in softening 
length which, although small, dominate over the errors due 
to timestepping. With the energy conserving term added 
(dashed line), only a small non-conservation of energy re- 
mains which can be shown to decrease as the timestep is 
made shorter. 

As a slightly more demanding test, we also consider the 
relaxation of a perturbed Plummer sphere - that is, with the 
velocities drawn from the equilibrium distribution function 
as described above, but then multiplied by a factor of 1.2. 
This means that the halo initially expands before slowly re- 
laxing into a dynamical equilibrium state. The evolution of 
the total energy in this case is shown in Figurc|S| Using adap- 
tive softening lengths, the change in softening lengths cor- 
responding to the initial expansion are reflected as a secular 
increase in the total energy (solid line) . Using the new energy 
conserving formalism (dashed line), this secular increase is 
not present and total energy is conserved to timestepping 
accuracy. 



4.4 SPH tests 

4-4-1 Static structure of a Polytrope 

A simple test of self-gravitating gas dynamics is to verify the 
static structure of a polytrope by allowing an initial arrange- 
ment of gas to settle into hydrostatic equilibrium. In order 
to do so we setup ~ 1000 SPH particles in a quasi-uniform 
spherical distribution and damp them into an equilibrium 
state using a polytropic equation of state P = Kp 1 with 
7 = 5/3. The low resolution is chosen in order to highlight 
the differences between various softening formalisms. 

The exact manner in which the particles are initially 
setup is not particularly important, although a perfectly 
uniform arrangement tends to produce numerical artifacts 
in the collapsed particle configuration whilst a clumpy ini- 
tial setup takes longer to settle to equilibrium. In this paper 
we use a quasi-uniform distribution achieved by placing par- 
ticles initially on a uniform square lattice, cropped to ensure 
that r < 1 and with a small, random perturbation of ampli- 
tude 0.2A (where A is the lattice spacing). The particle con- 
figuration is shifted slightly to ensure that the centre of mass 
is placed at the origin. Using a lattice spacing of A = 0.15 
results in a total of 1086 particles in the calculations. 




adoptive softening lengths (no energy conservation) 

adaptive softening lengths (with energy conservation) 
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Figure 5. Total energy conservation during the dynamical evo- 
lution of the 1000-particle Plummer sphere. Using adaptive soft- 
ening lengths without the additional term (solid line) leads to 
fluctuations in the total energy which dominate over the timestep- 
ping errors. Incorporating the new adaptive softening length term 
(dashed line), energy conservation is restored to timestepping ac- 
curacy. 




time 

Figure 6. Total energy conservation during the dynamical re- 
laxation of the perturbed 1000-particle Plummer sphere. In this 
case the initial velocities were multiplied by a factor of 1.2. Using 
adaptive softening lengths but without the new term (solid line) 
the change in the softening lengths caused by the initial expan- 
sion can be seen to cause a secular increase in the total energy. 
Adding this term (dashed line), the total energy is conserved to 
timestepping accuracy. 



The exact solution for the polytrope static structure is 
computed by solving the equation 



-yK 



d 2 



4ttG(7 - 1) dr 2 



[rp J 1 ] + rp = 0, 



(56) 



numerically using a simple finite difference scheme. The so- 
lution is then scaled to give a polytrope of radius unity. In 
code units (discussed in H4.3H a polytrope of radius unity is 
obtained by choosing K = 0.4246 in P = Kp 1 . 

In all simulations the density and SPH smoothing 
length are calculated by direct summation using the iter- 
ative method described in H4.2I Also the variable smoothing 
length terms in the SPH equations are used throughout. In 
order to isolate the effects of the softening formalisms we 
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Figure 7. Static structure of the 7 = 5/3 polytropc calculated 
using 1086 SPH particles (solid points). The results are shown 
using fixed plummer softening with softening length h = 0.03 
(top left), fixed cubic spline softening with h = 0.06 (top right), 
using adaptive softening lengths (bottom left) and finally using 
the energy-conserving formalism including the additional force 
term (bottom right). The exact solution is given by the solid line. 
Note that the SPH smoothing length is adaptive in all cases 



calculate the gravitational term by a direct summation over 
the particles (rather than using a treecode) with a standard 
second-order leapfrog scheme for time integration using a 
timcstep contr olled by a Coura nt condition based on the 
signal velocity jMonaghanll2005l) . The particles are damped 
to an equilibr ium using a stan dard form of the SPH artifi- 
cial viscosity ( Monaehan 119971) together with a damping in 
the force equation which is independent of particle number, 
given by 

^- = -0.05v + f , (57) 
at 

where f is the force per unit mass. Note that the polytropic 
equation of state means that the kinetic energy removed by 
the artificial viscosity and damping terms is not deposited 
as thermal energy but rather allowed to escape from the 
system. 

The equilibrium configuration of the 7 = 5/3 poly- 
trope with various softening formulations are shown in 
Figure [7] and may be compared in each case to the ex- 
act solution given by the solid line. The fractional errors 
(fi ~ f exact) / f exact are also shown in an inset plot in each 
panel. The top two panels show the results using fixed Plum- 
mer (top left) and cubic spline (top right) softening, where, 
not knowing the 'optim al' choice a prio ri we have used the 
rule-of-thumb given by ISpringe] i2005l) whereby the soft- 
ening length is chosen to be ~ 1/40 of the average parti- 
cle spacing in the initial conditions. Thus guided we choose 
h = 0.06 for the cubic spline softening, using half of this 
value, h = 0.03, in the Plummer softening (see discussion in 

urn . 

Using adaptive softening lengths without the energy 
conservation term (bottom left panel) shows a small im- 
provement over the fixed softening results, mainly in the 
outer regions where the force estimate is much less noisy. 
The density resolution in the centre is slightly lower in this 
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Figure 8. Total energy conservation during the radial oscillations 
of the polytrope. The results are shown using a fixed softening 
length (solid line), adaptive softening lengths (dashed line) and 
using the new adaptive softening length formalism (dot-dashed 
line). Note the improvement in energy conservation in the adap- 
tive softening case when the new term is included. The absolute 
value of the total energy differs slightly between runs because of 
the difference in equilibrium structure (Figure 0. 



case, but this is substantially improved when the energy 
conserving term is incorporated (lower right panel). The 
error in the outer regions is also improved by the energy- 
conservation term. The more compact distribution produced 
in this case is consistent with the additional term being al- 
ways in the direction of increasing the gravitational force 
(see ©. 

4-4-2 Polytrope oscillations 

Having obtained the static structure, studying the radial 
oscillations of the polytrope provides a test of the energy 
conservation properties of the code. In order to do so we ap- 
ply a radial velocity perturbation of v r = 0.2r to the static 
solutions obtained in the previous section. In order to dis- 
tinguish effects due to the softening formulation from effects 
due to the timestepping algorithm we use a very low Courant 
number of C CO ur = 0.05 for this test. In general, however, 
non- conservation effects from the softening formalism are 
much larger than effects due to timestepping. No artificial 
viscosity or damping is applied for this problem. 

The evolution of the total energy of the system is shown 
in Figure|H|using cubic spline softening with fixed and adap- 
tive softening lengths. Using a fixed softening length (solid 
line) the total energy is conserved exactly (ie. to timestep- 
ping accuracy). Adapting the softening length using the 
method described in H4. 21 results in non-conservation of en- 
ergy (dashed line). Incorporating the additional pseudo- 
pressure term into the adaptive softening formulation re- 
stores the total energy conservation (dot-dashed line). 



5 SUMMARY 

In this paper we have described an algorithm for using adap- 
tive softening lengths in both SPH and N-body codes which 
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retains the conservation of both momentum and energy. The 
formalism requires the computation of an additional grav- 
itational force term which is similar in form to the SPH 
pressure force and is therefore straightforward to implement 
in any SPH code at almost no added cost. For pure N-body 
codes calculation of the additional term requires some extra 
work since quantities such as the density must be evaluated 
using the smoothing kernel. However even in this case the 
cost is small compared to the evaluation of the long-range 
gravitational forces using a treecode. 

The softened gravitational force can be symmetrised ei- 
ther by using an average of the softening lengths or alterna- 
tively an average of the softening kernels, where the latter 
is preferred because of the manner in which the density is 
calculated. The choice of softening kernel is completely ar- 
bitrary, with calculations in this paper were made using the 
standard SPH c ubic spline kern el (although any of the ker- 
nels proposed bv lDehnenll200ll could be used). 

Use of spatially variable ('adaptive') softening lengths is 
found to provide near-optimal softening for arbitrary mass 
distributions using a single, fiducial choice of the adaptive 
softening parameter N ne j a h- This contrasts to the results of 
lAthanassoula et alJ <|2000T) where the optimal (fixed) soften- 
ing length was found to depend strongly on the number of 
particles and parameters such as the central concentration 
and shape of the mass distribution. For a mass distribution 
where more than one component is present, we find that the 
use of our adaptive softening length formalism can give more 
than an order of magnitude improvement in the errors on 
the least dense component. 

The main advantage of the formalism presented here 
is that adaptive softening lengths can be used whilst main- 
taining energy conservation to timestepping accuracy. This 
was found to be particularly important in the case of colli- 
sionless TV— body simulations where secular increases in the 
total energy were found to result from the use of adaptive 
softening lengths without the energy-conserving term. For 
self-gravitating SPH simulations the new formalism is a nat- 
ural and self-consistent choice which is found to give a small 
improvement in resolution and energy conservation over tra- 
ditional ad-hoc formulations for essentially zero additional 
cost. 
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APPENDIX A: CUBIC SPLINE SOFTENING 



In this appendix we give the functional form of the softening corresponding to the cubic spline kernel Hill . Integrating the 
kernel according to l|12[l . we find that the gravitational force is softened using 



<j>\r,h) = 

where q - 
given by 

4>(r,h) = 



(Al) 



1/h 2 , 0<g< 1; 

l/^ffg-S^ + f? 3 -^ 4 -^], 1< 9 <2; 
1/r 2 q > 2. 

r/h. Integrating a second time using 1131 gives the kernel used in the gravitational potential, which in this case is 
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< q < 1; 

Q > 2. 



The derivative of the potential with respect to h is given by 



dh 
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1/h 2 
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< q < 1; 

1 < g < 2; 
<7 > 2. 



Alternatively d<f>/dh can be evaluated from the potential and force functions according to 
where = /i</> and K' (q) = /i 2 <// are the functional forms of the potential and force kernels. 
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APPENDIX B: ADAPTIVE SOFTENING LENGTH FORMALISM USING AVERAGED SOFTENING 
LENGTHS 

A alternative way of symmetrising the gravitational potential is to use an average of the particle softe ning lengths. This is 
similar to the approach taken in the adaptive smoothing/softening length formalism used bv lBena (|199fl) . The main difference 
is that we symmetrise the gravitational potential rather than the force and are therefore able to account for the spatial variation 
of softening length in the equations of motion, leading to the conservation of both momentum and energy. Note that this is 
only made possible because of the self-consistent relationship between the density and the smoothing length described in H4.2I 
Using an average of the softening lengths, the gravitational part of the Lagrangian can be written in the form 

-— ^^m(,m c i , c (/i f , c ) (Bl) 

b c 

where (f> bc refers to (f>(\r b — r c |) and h bc = \(h b + h c ). It is then a straightforward matter to derive the equations of motion by 
using IB II in the Euler-Lagrange equations (1151 . The derivative of 1B1I is given by 
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Using the spatial derivative of the density given by 1231 and a similar expression for dp c /dr a we have 
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Collecting terms and simplifying, this expression can be written in the form 
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giving the TV— body equations of motion in the form 
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where in this case we define the quantity ( according to 
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This term is again easily calculated alongside p and Q during the density summation. However this formalism is quite 
inefficient in general, since the average h is only used in the calculation of The density, SPH and gravity forces are naturally 
symmetrised by the formulation from a Lagrangian. Calculation of £ in this case would require an extra loop over the particles. 
This is for the reason that, whilst the density and smoothing length can be iteratively found for a single particle a (depending 
only on h a ), quantities depending on an average smoothing length must be updated when a neighbouring value of h has 
changed (ie. ht), leading to a rather inefficient scheme . 

It is worth noting that lHernauist fc Barnes! il990f) suggested using a Lagrangian to derive an energy-conserving adaptive 
softening length formalism using an average of the softening lengths some time ago. However contrary to their assertion that 
'the terms involving Ve will, in general, lead to a violation of linear and angular momentum conservation', the force expressed 
by 1B5II clearly conserves linear momentum as the summations are antisymmetric in a and b. It is also straightforward to 
show that angular momentum is conserved exactly. 



